EFFECTIVE  METHODS  FOR  SOLUTION  OF 
NONLINEAR  REACTOR  DYNAMICS 
PROBLEMS  USING  FINITE 
ELEMENTS 


Richard  A1 len  Olsen 


asKssaSsr 


1 


JH 


ORterey,  uaiitorma 


^WjeJ^     B     ^M?^- 


Effective  Methods  for  Solution  of  Nonlinear  Reactor 
Dynamics  Problems  Using  Finite  Elements 


by 


Richard  Allen  Oisen 
December  1975 


Thesis  Advisor: 
Co-Advisor; 


D.H,  Nguyen 
D.  Salinas 


g^7»HIWll    jtt.i,  ».;* 


Approved  for  public  release;  distribution  unlimited. 


T170813 


UNCLASSIFIED 


SECURITY   CLASSIFICATION   OF    THIS  PAGE  flThan  Data  Entered) 


JWVAL  POSTGRA 
"ONTEREV,  CAU 


REPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


TTrepoRT  NUMBER 


2.  GOVT   ACCESSION  NO 


3.     RECIPIENT'S  CATALOG  NUMBER 


4.     TITLE  (and  Submit) 

Effective  Methods  for  Solution  of  Nonlinear 
Reactor  Dynamics  Problems  Using  Finite 
Elements 


S.     TYPE  OF   REPORT  b   PERIOD  COVERED 

Masters'  Thesis;  December 
1975 


6.  PERFORMING  ORG.  REPORT  NUMBER 


7.  AUThOR(«J 

Richard  Allen  Olsen 


6.     CONTRACT  OR  GRANT  NUMBEFlfaJ 


9.     PERFORMING  ORGANIZATION   NAME  AND  ADDRESS 


10.     PROGRAM  ELEUF.NT,  PROJECT,  TASK 
AREA  A  WORK  UNIT  NUMBERS 


Naval  Postgraduate  School 
Monterey,  CA   93940 


II.     CONTROLLING  OFFICE  NAME   AND  ADDRESS 

Naval  Postgraduate  School 
Monterey,  CA   93940 


12.  REPORT  DATE 

December  1975 


IS.  NUMBER  OF  PAGES 


(4.     MONITORING  AGENCY  NAME  A    ADDSSESSff/  dl  iterant  from  Controlling  Otllce) 

Naval  Postgraduate  School 
Monterey,  CA   93940 


IS.     SECURITY  CLASS,  (ol  thle  riport) 

Unclassified 


IS«.     DECLASSIFICATION/ DOWN  GRADING 
SCHEDULE 


16.     DISTRIBUTION   STATEMENT  (ol  thti  Report) 


Approved  for  public  release;  distribution  unlimited 


17.     DISTRIBUTION  STATEMENT  (ol  the  abrtrmct  entervd  In  Block  20,  II  different  from  Htpcri) 


ie.    supplementary  notes 


19.     KEY  WORDS  (Continue  on  revere*  aide  II  naceaecr?  end  Icenttty  by  block  nusnber) 


20.     ABSTRACT  (Continue  on  reverie  eldv  It  neceeaary  and  identity  by   block  nwrzber) 

The  solution  of  the  nonlinear  two-dimensional  reactor  dynamics  equation 

subjected  to  prompt  feedback  conditions  using  the  finite  element  technique 

leads  to  the  matrix  formulation  A.. v.  =  B  .  iK  +  C  .,  iK'Jj,  (iii.k  =1,  ....  N) . 

i]  J    ij  J    ijkr;Tk  '-" 

This  system  has  been  solved  directly  in  a  previous  work;  but  because  the 
nonlinearity  c--;t,^-j\  was  premultiplied  by   [A]-1,  large  computer  storage 


DD     t  JAN^J     1473  EDITION  OF    1  NOV  6»  IS  OBSOLETE 

(Page      1)  S/N    0102-314-6601   ! 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Vha,-.  Data  Cnterodi 


UNCLASSIFIED 


JtCUWTY    CLASSIFICATION    OF    THIS    PAGEfWi«n    P 


flm    £nlnf4i 


was  required  for  the  small  problem  considered.   The  task  of  this  thesis 
is  the  development  of  computational  techniques  which  allow  the  problem 
to  be  solved  for  large  systems.   Specifically,  these  techniques  are: 
(1)  the  treatment  of  the  nonlinearity  on  the  element  level,  (2)  the 
compacting  of  the  sparce  matrices  to  include  only  non-zero  terms,  and 
(3)  the  construction  of  a  new  computer  code  based  on  the  Crank-Nicolson 
formulation  for  the  solution  of  differential  equations. 

To  support  the  theory  presented,  test  problems  were  solved  by  the 
original  method,  the  linearized  technique,  and  the  Crank-Nicolson 
treatment.   The  results  were  analyzed  and  compared  graphically.   All 
three  of  the  innovations  developed  in  this  thesis  appear  to  be  useful 
tools  for  solving  nonlinear  time  dependent  differential  equations. 


^lKmSa  U73  UNCLASSIFIED 

S/  N       0102-014-6601  SECURITY   CLASSIFICATION   OF   THIS   PAGEr***"  Dmf   Bn*er.d) 


\ 


Effective  Methods  for  Solution  of  Nonlinear  Reactor 
Dynamics  Problems  Using  Finite  Elements 


by 


Richard  Allen  Olsen 
Lieutenant  Commander,  United  States  Navy 
B.S.,  United  States  Naval  Academy,  1966 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  MECHANICAL  ENGINEERING 

from  the 

NAVAL  POSTGRADUATE  SCHOOL 
December  1975 


OUDLEY  KNOX 

NAVAL  POSTGRADUATE  SCHOOL 

MONTEREY.  CALIFORNIA  93940 


ABSTRACT 

The  solution  of  the  nonlinear  two-dimensional  reactor  dynamics 

equation  subjected  to  prompt  feedback  conditions  using  the  finite 

element  technique  leads  to  the  matrix  formulation  A .  .  iL  .  =  B..iK  + 

ij  J    ij  j 

C..,i^.^,  (i,j,k  =  1,  ...»  N) .   This  system  has  been  solved  directly 
in  a  previous  work;  but   because   the  nonlinearity   C.,,  ty.ty,      was 
premultiplied  by   [A]   ,  large  computer  storage  was  required  for  the 
small  problem  considered.   The  task  of  this  thesis  is  the  development 
of  computational  techniques  which  allow  the  problem  to  be  solved  for 
large  systems.   Specifically,  these  techniques  are:   (1)  the  treatment 
of  the  nonlinearity  on  the  element  level,  (2)  the  compacting  of  the 
sparce  matrices  to  include  only  non-zero  terms,  and  (3)  the  construction 
of  a  new  computer  code  based  on  the  Crank-Nicolson  formulation  for  the 
solution  of  differential  equations. 

To  support  the  theory  presented,  test  problems  were  solved  by  the 
original  method,  the  linearized  technique,  and  the  Crank-Nicolson 
treatment.   The  results  were  analyzed  and  compared  graphically.   All 
three  of  the  innovations  developed  in  this  thesis  appear  to  be  useful 
tools  for  solving  nonlinear  time  dependent  differential  equations. 
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I.   INTRODUCTION 

When  temperature-dependent  feedback  is  considered  in  the  nuclear 

reactor  dynamics  problem,  a  nonlinear  field  equation  in  space  and  time 

results  [Ref.  1].   For  the  non-homogeneous,  or  multi-region  reactor, 

however,  the  space-dependent  dynamics  behavior  following  a  nonlinear 

initial  disturbance  is  no  longer  reachable  in  analytical  form. 

Fortunately,  by  modeling  the  reactor  as  a  system  of  finite  regions  of 

interest  where  the  neutronic  properties  are  known,  the  method  of  finite 

elements  can  be  applied  to  yield  the  solutions.   The  fundamental  concepts 

relating  reactor  behavior  to  the  finite  element  formulation  have  been 

presented  in  19  74  [Ref.  2],  but  a  more  recent  work  by  Nguyen  and 

Salinas  [Ref.  3]  gives  a  thorough  discussion  of  the  complete  problem. 

That  work  forms  the  basis  and  starting  point  for  this  thesis,  the 

objective  of  which  is  to  develop  improved  computational  methods  for 

dealing  with  that  finite  element  formulation. 

In  Reference  3  the  dynamic  flux  equation  under  prompt  feedback 
conditions  was  given  as 

.      3tp    (r,z,t)  .  9 

i    — 2- =  D  V  if)     +  y   Z      4     -   a  K  I      tf>  (1) 

V  3t  mm  mamm  mmamm 

m 

for  each  non-homogeneous  zone   "m"   contained  in  the  reacLor  body, 
where  the  usual  symbols  are  employed: 


X     =  v  l£   /I       -  1 
m      fm  am 


K  =  e  l£   /p  C    (°C/unit  flux-sec) 
m      fm  m  pm 


Ym  =  (Am/V  )  (hm/pmC   )    sec  X 
m     m  m   m  m  pm 


p  C  =  heat  capacity 

a  =  reactivity  temperature  coefficient  °C 

h  =  convection  heat  transfer  coefficient 

A/V  =  heat  transfer  area  to  volume  of  energy 
generation  ratio 

Ef  =  neutron  fission  cross-section 

£  =  neutron  absorption  cross-section 
a  r 

v  =  neutrons  emitted  per  fission 

e  =  energy  produced  per  fission 
The  subscript  "m"  will  be  omitted  from  further  discussion  for  simplicity. 
The  transformation  of  this  general  equation  into  a  problem  in 

finite  dimensional  vector  space  begins  by  constructing  the  following  N 

term  approximation 

N 


ij>(xft)  =  iC(x,t)   =  ?      <K(t)G,  (x) 


(2) 


where  N  is  the  number  of  degrees  of  freedom,  or  nodes,  in  the  space, 
and  the  G.  are  the  basis  functions  of  the  approximate  solution  space 
The  Galerkin  method  seeks  to  make  the  residual 

R(x,t)  ■  L  j  -  f , 

where  Liji  =  f   is  the  field  equation,  othogonal  to  each  of  the  basis 
functions  so  that 


This  discussion  leading  to  the  establishment  of  the  matrix 
equations  is  essentially  an  abridgement  of  the  development  given  in 
reference  3. 


/  G  (x)R(x,t)dx  =0    i  =  1,  2,  ...  ,  N  . 


(3) 


For  the  nuclear  reactor  dynamics  problem  of  Eq.  (1),  the  residual 
becomes 

2 

R(x,t)   -  !?■  -  V  DV2^  -  VAI  J  +  wE  I2 
—        3t  a      a 


(4) 


with  VaKE  =  w  .   Putting  this  expression  into  Eq.  (3)  and  integrating 
a 

by  parts   (a  distinct  advantage  of  Galerkin) ,  the  following  coefficients 
emerge 


IJ 


"/A 


Gv   rdrdz 

K. 


(5a) 


IJ 


S  L 


3G,    3GT      3GT    3G, 


+ 


3r    3r      3z    3z 


rdrdz 


(5b) 


CIJK  =  iI  GIGJGKrdrdz  (5C) 

I,J,K  =  1 N 


in  cylindrical  coordinates  where   dS  =  rdrdz.   This  allows  Eq.  (1)  to  be 
written  as  follows  for  the  case  of  uniform  neutronic  properties  within  each 


reactor  region  m. 

N 


N 


I  Au*j  ■ "  vi 


J=l 


D!   S  BIJ*J  +  VlEaI  Z  AIJ*J  " 

J=l  J=l 


N    N 
J=l   K=l 


wi  E  E  cijkVk    i  =  1>2>  '••  N 


(6) 
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This  becomes,  upon  combining  constants,  in  Einstein  summation  convention, 
where  there  is  no  sum  on  an  underlined  repeated  index, 

AIJ  *J  "  ABIJ  *J  "  WICIJK  *A 

(7a) 

*i  =  W^  I.JtKi  ■  li  •••»  N 

With  boundary  conditions,  the  system  is  now  well  posed  and  ready  for 
solution.   The  equation-solver  must  be  chosen  with  care,  however, 
because  this  system  of  ordinary  differential  equations  is  both  stiff 
and  non-linear.   "Stiff"  is  used  here  and  throughout  this  thesis  to 
denote  a  system  which  gives  a  large  response  to  a  small  stimulus;  in 

this  case,  the  nature  of  the  reactor  dynamics  problem  predicts  a  large 

2 
change  in  flux  over  a  small  change  in  time.    Current  experience 

suggests  Gear's  predictor-corrector  method  [Ref.  4]  as  written  in  the 

computer  programming  code  DVOGER  [Ref.  5].   It  requires  only   1)  the 

locus  (value)  of  the  points  at  a  given  time,   2)  a  routine  to  evaluate 

the  instantaneous  derivatives  at  any  time  and   3)  a  routine  to  evaluate 

the  Jacobian   (J  =  3f  .  (ij>,  t)/3^  .  )  .   Based  on  this  information,  the  program 

assumes  a  time  step,  which  may  be  quite  small  (the  minimum  size  being  a 

function  of  the  particular  computer) ,  and  tries  to  fit  a  trial  (predicted) 

solution  for  that  time  interval.   The  corrected  (integrated)  solution  is 

then  obtained  by  iteration  until  convergence  is  attained.   If  the 

predicted  and  corrected  solutions  fail  to  match  within  a  specified  error 

criterion,  the  time  step  is  decreased  and  the  process  repeated.   On  the 

other  hand,  if  "excessive"  accuracy  is  attained,  the  next  attempted 

time  step  will  be  increased.   Gear's  method  in  DVOGER  is  well  suited 

to  stiff  non-linear  systems.   Also  contained  in  DVOGER  is  the  Adam's 


2 

This  definition  is  somewhat  broader  in  scope  than  that  used  by 

many  other  authors. 
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method  which  parallels  that  of  Gear  except  that  the  Jacobian  is  not 
computed.   Lacking  the  additional  information  about  the  rate  of  change, 
the  Adam's  method  must  proceed  much  more  cautiously  with  stiff  systems 
and  hence  marches  forward  with  exceedingly  small  time  steps.   Present 
experience  confirms  the  caveat  advice  given  by  DVOGER  that  the  Adam's 


me 


thod  is  not  intended  for  "stiff"  cases. 


The  difficulty  that  arises  from  acquiring  solutions  to  Eq.  (7), 
then,  should  not,  and,  in  light  of  the  present  experience,  does  not 
result  from  DVOGER  or  any  other  acceptable  equation  solver,  but  is,  rather, 
a  function  of  matrix  size,  vis  a  vis  the  number  of  nodal  points  in  the 
finite  element  approximation.   The  finite  element  modeling  of  large 
reactors  or  the  close  scrutiny  of  small  ones  is  severely  restricted, 
therefore,  unless  the  number  of  nodal  points  can  be  raised  significantly. 
As  an  example,  the  test  case  considered  in  this  thesis  consisted  of  only 
38  nodes  but,  when  processed  directly  from  Eq.  (7),  required  over  300K 
bytes  of  storage  in  the  IBM  360/67  computer.   This  requirement  came 
largely  from  the  [A]   [C    ]  matrix  alone  needing  a  size  of   4(38x38x38)  = 
219K  bytes  in  single  precision.   Manipulations  with  this  cubic  were 
indeed  quite  burdensome.   Since  machine  processing  time  is,  other  things 
being  equal,  dependent  on  size,  the  combined  time  and  space  requirements 
would  prohibit  the  consideration  of  even  moderate  sized  problems  except 
on  the  largest  computers.   Additionally,  the  inversion  of   [A]   in  Eq. 
(7)  is  indicated  in  order  to  provide  an  explicit  value  of  \\>     for  the 
DVOGER  routine.   It  is  this  difficulty  of  size,  and  with  it  processing 
time,  to  which  the  remainder  of  this  thesis  is  devoted.   Experience  with 
the  DVOGER  routine  is  presented,  and  the  author's   equation-solver  based  on 


12 


a  Crank-Nicolson  formulation  is  discussed  as  the  analysis  of  the  size 
problem  is  developed.   Table  I  shows  clearly  the  influence  of  matrix 
size  on  computer  storage  requirements. 


13 


II.   REDUCTION  OF  THE  NONLINEARITY 

A.  GLOBAL  TREATMENT  OF  THE  NONLINEARITY 

The  evaluation  of  Eq.  (7)  generates  the  following  global  matrix 
system,  where  the  subscripts  indicate  the  matrix  size  required  for  N 
number  of  nodal  points 

[^xn]{*Nxi}   =   [BNxn]{W}  +   [CNxNxN]{'|,(NxN)xl}  (8) 

2 
where   i|)   is  the  time  response  of  the  flux  at  each  nodal  point  and  \\> 

is  the  product   ^.^.   i»j  =  1>2,  ...  N.   Solved  in  this  form,  the  major 

3 

limitation  of  the  procedure  was  the  large  N   size  of  the  cubic  array 

C,  which  is  introduced  as  the  result  of  the  nonlinearity  of  Eq.  (1). 

With  this  size  requirement  rapidly  increasing  for  a  finer  mesh  consisting 

3 

of  more  nodal  points,   the  investigation  of  Reference  3  was  limited  to  a 

small  but  practical  N   of  38. 

B.  TREATMENT  ON  THE  ELEMENT  LEVEL 

The  cubic  form  of   C  need  not  appear  if  it  can  be  shown  that  the 
nonlinear  portion  of  Eq.  (8)  can  be  broken  apart  into  the  product  of 
two  linear  components  such  that  one  component  can  be  evaluated  on  the 
element  level  and  the  element  contributions  then  summed  into  global  form, 


3 

In  general,  for  a  nonlinearity  of  order  m,  a  finite  element 

problem  of  N  nodes  will  generate  a  matirx  size  of  Nm  ^. 
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resulting  in  a  "regular"  NxN  matrix  C*  ready  for  multiplication  by 
the  remaining  ty  also  summed  to  the  global  level.  Proof  that  this  is 
indeed  the  case,  i.e., 


[CNx(NxN)|{,|;(NxN)xl}     [SxnJ  {^Nxl} 


(9) 


is  outlined  as  follows: 

Reference  3  established  the  element  nonlinear  term  as 


/A^vw^^ijkV0  *k(t) 


(10) 


where  the  integral  is  over  the  element  area  A   , 


Cijk  =  Jl     ^jCk(re^e)drdz'    ijj'k'e  =  1'2'3 


(11) 


and   £  being  the  local  triangular  coordinates,  an  innovation  of  Felippa 
[Ref.  6].   The  coefficients   C. .,   form  a  3x3x3  element  array  possessing 
a  regularity  which  suggests  a  kind  of  "cubic  symmetry."  The  evaluation 
of  this  integral  results  in  the  following  expression 


Clll  '  *"[*rl  +  6(r2  +  'P] 

re[eri 


112 

Cn3 =  Y 

C  =   y 

L122        Y 

C  -  y 

L123        Y 


+  Ar2  +  3r3] 


T6ri  +  2r2  +  4r3] 


C133  "  Y 


l[2(rx 

T4ri 


+  6r2  +  2r3] 

+  r2  +  r3>] 
+  2r2  +  6r3] 


'222 


Y6^^  +  24r2  +  6rJ 


C333        Y 


K 


+  6r2  +  24r3] 


(12) 
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where 


and 


Y  -  frA  /180 
e 


c         =  r    =  c 
121    112    211 


C321   C231  "  C213  etc  as  in  Ref.  3. 


2  2 

Assuming   that      ty       may  be  broken  apart   as      ty     =   \\>    •    <j>,    the  LHS 


(left  hand  side)  of  Eq.  (10)  is  rewritten  as 


ffwsw 


i,j,k  =   1,2,3 


(13) 


or 


27T         ff^     [?1C2?3] 


♦l 

*2 

*3 

[cxc2c3] 


♦l 

*2 

*3 

r  dr  dz 


(14) 


Expanding  and  collecting  terms  yields 


2tt 


//«k 


*1  *l(rlCl  +  W2  +  VlC3> 


+  ♦jC^cJtj  +  r2c1?2  +  r3?lC2C3) 


(15) 


+  +3(^3  +  r2ClC2C3  +  r3Cl^) 


+  *2  K(W2  +  r2^2  +  W2C3>   +  ^2(rri^2  +  rl4  +  r3^3} 
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+  VW2?3  +  W3  +  W3} 


+  ^   ^(r^^  +  V1?2C3  +  ^.C*)  +  ^(r^^  +  t^   +  ^C*) 


+  ♦3(r1C1C3  +  r^C*  +  ^ 


dr  dz 


k  =  1,  2,  3 


Since  the  quantities   if»   and   <J> .   (i»j  =  1,2,3)   depend  only  on  time, 
they  are  unaffected  by  the  integration.   Further,  the  expressions 
inside  the  parentheses  are  given  a  unique  name  d.  ,,   to  fix  their 
position  in  relation  to  a  particular  ty.fy.      f°r  some  k  =  1,2,3.   Now 
the  above  expression  may  be  written  as 


2tt  if>.d>.   //  d.  ...  dr  dz  . 
riYj  J  J       xjk 


(16) 


Using  the  integration  formula  from  Felippa  [Ref.  6] 


// 


I     m  n   ,        ft!m!n! 
Cl  C2  C3  dA  "   (Jl-hn+n+2)! 


x  2A 


(17) 


where   £,m,n  are  the  exponents  of   £   in  a  specific  d...  .   The  integration 

ijk 

of    (16)    gives 


(ttA  /180)    iK4>.    f .  .. 
e  iYj      ijk 


(18) 


and      f     ,    x  ttA  /180      is    identically   equal  to      C.  .,       of    (12). 
x  j  tc  e  13  *£ 

The   indicated  summation  is    conveniently  expressed  by 


<K    f .  ,,    -  a.        i,j,k  =   1,2,3 
j      ijk         x  J 


(19) 
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e       * 

which  is  immediately  recognized  as  y   a   =  C . ,  ,   the  desired  3x3 

* 

element  matrix.   In  this  form  it  is  combined  into  the  global   [C„  „]  , 

NxN 

and  the  linearization  has  been  accomplished.   As  will  be  shown  later, 

it  may  be  convenient  to  let   d>  =  \b      ,        when  for  small  At.'  6  B 

t-At 

1.   s  ^   .   This  approximation  is  useful  for  predictor-corrector 

equation  solvers  in  order  to  allow  construction  of   [C„  „1   only  once 

NxN 

for  each  successful  time  step. 
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III.   SOLUTION  PROCEDURE  WITH  THE  LINEARIZED  FORM 


For  clarity,  throughout  the  remainder  of  this  thesis  the  following 

3 

naming  convention  will  be  used:   N   denotes  the  nonlinearized  direct 

treatment  by  DVOGER  as  solved  in  Ref.  3.   "Linearized  approximate" 
indicates  the  linearized  technique  based  on  element  level  multiplication 
of   C  by  lPt._At'      The  "linearized  exact"  method,  however,  is  the  element 
level  multiplication  of   C  by  the  trial  ip   .   The  compact,  or  Nxp 
formulation  may  be  used  with  any  of  the  linearized  methods  in  addition 
to  the  CRANKO  treatment  discussed  later. 


A.    GENERAL  TREATMENT  OF  THE  MATRIX  EQUATION 

Solution  of  the  matrix  formulation  of  Eq.  (8)  involves  clearing 
the  LHS  to  provide  a  discrete  relation  for  each  ip,   .   This  is  accomplished 
either  by  iteration,  which  yields  an  approximation,  or  by  matrix  inversion, 

which  has  the  advantage  of  giving  the  individual  ip,       explicitly.   In 

3 
both  the  nonlinearized  N   and  the  linearized  NxN   forms,  Eq.  (8)  is 

reduced  with  multiplication  by   [A]   ;  thus 

[AfVui)   =    [A]_1[B]{^}  +    [A]_1[C<j)]{^}  (20) 

2 
in   the   linearized      \p     =   $\p      case   or   as 

[AfVui}   -    [A]"1[B]{lf)}   +    [Af^CJU2}  (21) 


for   the   direct,    nonlinear,    formulation.      This   results   in 


W   =    [AB]{^}  +    [A   1][C*]{^}  (22) 
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where 

[C*]  =  [C4>]   or  [C)W  (23) 

as  appropriate.   The  point  here  is  that  the   [AB]   and  the   [A   ]   are 
formed  once  and  for  all  based  on  the  time  independent  properties  and 
geometry  of  the  problem,  whereas  the  [C$]      must  be  computed  for  each  time 

t.   Although  the  construction  of  the   [C]   is  required  only  once  in  the 

3  2 

direct   (N  )   method,  it  must  be  additionally  multiplied  by  {$    }      at 

each  time  step. 


B.    PRACTICAL  CONSIDERATIONS  WITH  THE  LINEARIZED  TECHNIQUE 

Relating  Eq.  (22)  to  the  demands  of  the  subroutine  DVOGER,  which 
involves  the  evaluation  of   ip   for  every  trial  value  of   i/j,   it  is 
apparent  that  the  construction  of   [A]   [C^]{i^}   for  each  of  these  trials 
could  be  a  time  consuming  process,  and  hence  a  serious  drawback  of  the 
linearized  treatment.   Noting,  however,  that  for  small  At,$  s&  i/>.  Al_> 
the  approximation   [Ci>]  ^  [Cip ]  ,   where   <f>  =  i);    ,  overcomes  this 
disadvantage  by  allowing  the  construction  of   [A]   [C<f>]   only  once  for  a 
given  time   t,  regardless  of  how  many  trial  solutions  are  attempted  by 
DVOGER.   Equation  (22)  can  now  be  further  simplified  since   [A]   [Ct}>] 
is  readily  combined  with   [AB]   to  reformulate  the  problem  as 

W   =    [C**]{*}  .  (24) 

Thus,  although  the  linearized  form  can  be  handled  exactly,  i.e.,  with 
4>  =  i|>,  it  may  be  more  expeditous  to  establish  the   [C   ]   only  once  per 
valid  time  point,  (that  is,  where  the  convergence  criterion  of  the  equation 
solution  has  been  satisfied).   In  this  case,  the  linearized  treatment  is 
then  an  approximate  technique  in  that  $   =  i(i   .  4   tyt> 
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C.    DIFFICULTIES  WITH  THE  APPROXIMATION 

The  simplification  expressed  as  Eq.  (24)  works  extremely  well  when 

used  with  DVOGER  as  long  as  large  changes  in  ib     are  incurred  by  small 

At   (the  so-called  "stiff"  system).   When  the  solution  approaches  steady 

state,  however,  ip      approaches  zero  faster  than  (ip     -  4>  )/At   does, 

and  this  causes  a  second  kind  of  perturbation  wherein  the  system  becomes 

increasingly  sensitive  to  small  changes  in  ip     -   ip   _.. •   Apparently 

DVOGER  requires  very  accurate  ip     information.   The   tj»(t)   provided  from 

i>  /    \  x  ip ,      ._*   does  not  match  the  DVOGER  prediction  of  ip,J_si   which  is 
r(t)     (t-At)  r  (t) 

2 
based  on  the  assumption  that   if),  .   was  generated  by   ^/^s*   This 

discrepancy  in  the  neighborhood  of  ib   *  0   causes  DVOGER  to  be  unable  to 

recognize  the  steady  state  solution.   Prepared  especially  for  stiff 

systems,  the  equation-solver  expects,  instead,  a  large  scale  change  in 

ib.      The  time  step  is  therefore  reduced  accordingly,  which  means  that  very 

slow  progress  is  made  in  the  steady-state  region  which  is  otherwise 

handled  quite  rapidly  under  the  exact  methods.   The  implication  here  is 

that,  in  those  cases  which  reach  terminal  flux  values  rather  early  in 

problem  time,  the  exact  formulation  of  the  linearized  treatment  is  to  be 

preferred,  even  though  the   [C   ]   must  be  computed  for  each  attempted 

ib.   ..   This  aspect  has  been  analyzed  for  the  test  problems  considered, 

with  the  results  given  in  Table  II. 
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IV.   THE  Nxp  COMPACTING  SCHEME 

A.    ESTABLISHED  COMPACTING  METHOD 

The  problem  under  consideration  here,  being  of  the  form  A.  ,i|>,  = 
B.(^',  ...  ifO ,  is  formally  handled  by  multiplying  through  by   [A]   , 
discussed  earlier.   For  large  matrices,  the  inversion  process,  and  indeed 
all  matrix  operations,  become  extremely  costly  in  terms  of  machine  time 
and  storage.   The  fact  that  the  matrices  may  be  banded  or  symmetric  does 
offer  significant  economy  when  the  matrix  is  not  inverted.   For  example, 
a  symmetric  NxN  matrix  may  be  stored  as  NxN/2.   For  matrices  resulting 
from  finite  element  formulation,  the  numbering  schedule  determines  a  band 
width  to  be  used  with  conventional  compacting  schemes.   In  the  finite 
element  case,  the  difference  is  compared  for  each  node  in  the  mesh.   The 
largest  of  these  differences  plus  one  equals  the  minimum  band  width 
allowed.   For  example,  node  7  in  Figure  1  has  a  maximum  difference  of 
13  -  2  =  11,  whereas  node  22  has  the  largest  difference  (37  -  10  =  27) 
for  the  system.   Thus  the  resultant  minimum  band  width  would  be  28.   As 
the  mesh  increases  in  nodal  points,  the  minimum  band  width  must  also 
increase  in  response,  no  matter  how  adroitly  the  numbers  are  assigned. 
Regardless  of  how  sparce  a  matrix  may  be,  these  schemes  all  ensure  that 
the  "compactness"  grows  commensurately.  In  addition,  the  calculation  of 
an  inverse  becomes  ever  more  unpleasant,  especially  since  this  operation 
on  a  sparce  matrix  generally  results  in  a  dense,  nonsymmetric  inverse. 
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B.    BASIS  FOR  THE  COMPACTING  SCHEME 

In  solving  a  matrix  equation  A..il>.  =  B.(ii.  ,  ...  ,  ib  )   it  is  not 

ij  J     i  1         N 

necessary  to  form  the  inverse,  but  rather,  multiply  out  the  left  hand 
side  and  rearrange  to  obtain  a  set  of  equations: 


<h  =  (B.-A.  .ijO/A..  (25) 

i     i   ij  j    ii 


where   i,j,k  =1,  ...  N.   This  procedure  is  especially  suitable  for 

sparce  systems  since  the  number  of  non-zero  A. .   entities  will  be  small. 
'  ij 

Further,  if  these  non-zero  A.,   can  be  located  and  tagged  with  an 
identifier,  then  they  may  be  stored  together  in  a  dense  array,  thereby 
replacing  the  standard  NxN  size  matrix  by  an  equivalent  array  of  size 
Nxp,  where  p   is  the  maximum  number  of  non-zero  entries  on  any  one  row. 

This  type  of  compacting  is  especially  well  suited  to  matrices 
resulting  from  finite  element  formulation,  with  p   determined  by  the 
choice  of  the  finite  element  and  the  discretized  model.   In  the  case 
under  study,  a  linear  triangular  element  is  used  such  that,  when  assembled 
on  the  global  level,  it  forms,  for  the  purpose  of  illustration,  a  network 
of  interlocking  hexagonal  polygons,  each  having  an  apex  located  directly 
over  its  parent  central  node.   Such  a  pattern  is  sketched  as  Figure  2. 
Note  that  the  sides  are  indeed  the  contributing  "neighbor"  nodes.   Boundary 
nodes  will,  of  course,  be  missing  any  "would  be"  contributors  sought  in 
the  boundary  exterior.   In  a  regularly  drawn  system,  as  in  Figure  2,  an 
interior  node  will  have  six  neighbors  plus  itself  for  a  total  of  seven 
contributors.   That  is  to  say,  that  no  matter  how  large  the  system  of 
mesh  points,  a  relation  describing  the  effect  of  the  system  on  any  given 
point  would  contain  a  maximum  of  seven  non-zero  values. 
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Geometric  considerations  or  a  desire  to  examine  some  area  of  the 
mesh  more  closely  may  result  in  systems  of  irregular  discretization 
wherein  the  maximum  number  of  entries  may  be  greater  than  seven.   This 
analysis  may  be  extended  to  consider  other  finite  element  shapes  in  a 
similar  manner. 

C.    CONSTRUCTION  OF  THE  Nxp  ARRAY 

To  utilize  this  attribute  of  finite  elements  to  produce  the  dense 
Nxp  system,  it  is  necessary  only  to  construct  a  table  of  nodal  points 
listing  their  contributors.   This  table,  or  connectivity  array,  is  then 
used  as  an  index  to  locate  the  values  of  the  various  quantities  associa- 
ted with  a  particular  node.   All  the  standard  matrix  operations  can  be 
performed  on  these  compacted  arrays,  but  the  construction  of  a  matrix 
inverse  has  no  usefulness  since  a  NxN   system  is  then  reconstituted  from 
an  Nxp   array. 

As  an  example,  construction  of  the  connectivity  vector  associated 
with  node  18  in  Figure  1  is 

[18  17   27   26  19  14  13  12] 

and  the  vector  associated  with  node  6  is 

[6  11   7   21  0  0  0  0]. 

The  contributors  may  be  entered  in  any  order  except  that,  for  computa- 
tional ease,  the  central,  or  "parent"  node  is  the  first  element.   These 
vectors  are  collected  into  the  node  neighbor  connectivity  array  of  size 
Nxp,  or,  for  Figure  1,  38x8. 
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To  multiply  a  Nxp   array  of   a..   elements  by  a  vector  of  b. 
elements,  a  search  of  the  connectivity  array  is  performed  on  the  "i"th 
row  to  locate  a  match  for  the  b..   When  the  match  is  found,  the  values 
represented  by   a..   and  b.   are  multiplied  together  resulting  in  a  new 
C.   vector.   If  the  system  is  correctly  formed  and  compacted,  there 
always  will  be  a  match. 

The  compacting  scheme  is  bound  by  the  geometry  of  the  problem,  and 
as  such,  has  meaning  only  as  it  pertains  to  that  problem.   For  example, 
operations  with  two  equi-sized  Nxp   arrays  of  different  connectivity 
relations  cannot  be  performed.   This  is  in  contrast  to  regular  matrix 
operations  where  the  restriction  is  only  to  size  and  not  to  origin.   A 
short  computer  program  to  demonstrate  the  operation  of  this   Nxp   scheme 
is  given,  with  results,  in  Appendix  C. 
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V.   THE  CRANK-NICOLSON  FORMULATION 

Recognizing  the  very  sparce  nature  of  the  system  of  differential 
equations  formulated  by  this  problem  and,  if  the  matrix  inversion  process 
is  to  be  avoided,  the  necessity  of  iteration,  a  straightforward  attack 
based  on  the  definition  of  the  derivative  (finite  difference)  appears  as 
a  feasible  alternative  to  the  DVOGER  technique.   One  such  approach  was 
presented  by  J.  Crank  and  P.  Nicolson  in  1947  [Ref.  7]  and  has  been 
discussed  in  many  works  throughout  the  ensuing  years.   This  method  has 
been  shown  to  be  unconditionally  stable.   An  adaptation  of  this  method 
is  given  here  as  follows: 

In  the  general  sense,  a  matrix  formulation  of  a  set  of  linear 
differential  equations,  can  be  represented  as 


A.  .$.    =  C.  .i>.   +   F.  (t) 
13  J     i]  ]     J 


i,j  =  1,  . . .  ,  N 


(26) 


where   ^   is  a  function  of  time.   For  the  case  of  an  initial  disturbance 
as  the  forcing  function,   F(t)  =  0  for   t  >  0.   Focusing  attention  for 
the  moment  on  a  particular  equation  and  expressing  [16]  in  terms  of  the 
definition  of  the  derivative 


which  is 


yt  yt-At 
At 


=   C 


yt    rt-At 


(27) 


[A-   1/2'Af  C]i|/  =  [A  +  (1/2-Af  C)]*t_At  (28) 


[(2/AfA)-CHt  =  [(2/AfA)+C]^t_At 


(29) 
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letting   D  =  (2/AfA)-C,   E  =  (2/AfA)+C,   <j>  =  ^t_Afc   then,  returning 
to  the  system  form  the  result  is 

Dij^j  =  Eij*j  =  Ei*    i.J  =  1»  .-»  N  .        (30) 


This  constitutes  a  set  of  linear  simultaneous  equations  which  may  be 
solved  by  any  convenient  method.   Chosen  here  is  Gauss-Seidel  iteration 
[Ref.  8]  which  is  well  suited  for  the  problem  at.  hand  since,  the  sum  of 
D.  .iK   and  E..d>.   will  contain  very  few  non-zero  terms,  thereby  incurring 
relatively  small  roundoff  errors  regardless  of  system  size. 

This  procedure,  which  also  uses  the  Nxp   compacting,  has  been 
written  for  this  thesis  as  the  computer  code  CRANKO.   After  receiving 
the  matrix  information  and  control  parameters  from  the  calling  program, 
CRANKO  first  builds  the   C  matrix  of  Eq.  (26)  for  the  initial  <\>     value. 
Then  it  selects  a  trial  time  interval  and,  forming  the  equation  set  (30) 
based  on  this   At,  attempts  to  find  a  solution  satisfying  a  specified 
convergence  criterion  based  on  a  relative  error.   If  no  satisfactory 
solution  set  is  found  after  a  number  of  iterations,  a  smaller  At   is 
selected;  and,  after  recomputation  of  Eq.  (30),  another  attempt  at 
convergence  is  made.   The  successful  solutions  then  replace  the  previous 
ty     thereby  forming  the  starting  point  for  another  cycle,  beginning  with 
a  new  computation  of  the   C  matrix.   The  choice  of   At   is  the  control- 
ling factor.  The  scheme  employed  counts  the  number  of  iterations  required 
for  a  solution.   When  an  increasing  number  of  iterations  indicate  greater 
difficulty,  the  interval   At   is  decreased;  whereas,  if  fewer  iterations 
are  required,  the  time  step  may  be  relaxed  and  made  larger.   Experience 
with  the  CRANKO  routine  establishes  this  method  as  an  extremely  fast 
technique  for  solving  this  type  of  problem. 
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For  those  readers  unfamiliar  with  the  Gauss-Seidel  iteration  scheme, 
a  simple  example  using  a  3x3  system  is  given: 


D  .$.    -  EJ     6   =    ih 


"dll   dl2   d13" 


d21   d22   d23 


d31   d32   d33J 


r*i_ 

~er 

*2 

= 

e2 

u_ 

_e3_ 

first  iteration: 


dll^l  +  d12*2  +  d13*3  =  el 


=    [e,    " 


=   *i 


(diA  +  diA)]/dn 


d21*1  +  d22^  +  d23^ 


=   e, 


=    [e0  " 


=   *. 


(d2A  +  d23*3)]/d22 


d31^1  +  d32*2  +  d33*3 


=  e. 


*3=  te3-  (d3l*l+d32*2)!/d33 


for  successive  iterations, 
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dll*l  +  d12*2  +  d13*3  =  el 


♦l"    [el"    (d12*2+d13*3>1/dll 


d21*l  +  d22*2  +  d13*3  =   e2 


•k-k  A 

*2  =    [e2  "    (d2lA     +  d13*3)]/d22 


2  ''       2 


** 


d31*l  +  d32^2  +  d33*3  =  e3 


etc. 


♦3-   [e3-   (d31*1    +d32*2)]/d33 
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VI.   SUMMARY  OF  AVAILABLE  METHODS 

For  clarity,  the  innovative  computational  tools  introduced  in  this 
thesis  are  summarized.   First  is  the  linearizing  of  the  governing  equation 
Eq.  (8)  by  the  multiplication  on  the  element  level  of   C...   by  \p   . 
From  this  development,  two  options  appear  feasible:   either  the  exact 
computation  of   [C  ]   for  each  trial  ty,    or  the  approximation  of 
\l>     st  ^r_Ar      to  allow  the  construction  of   [C  ]   only  once  during  the 
search  for  a  new  ty   .      This  approach  uses,  as  does  the  "direct,"  or 
nonlinear  treatment  involving   [C],  the  inverting  of   [A]   to  clear  the 
LHS  as  shown  in  Eq.  (22). 

The  second  innovation  is  the  reduction  of  the  NxN  matrix  to  a 
smaller  Nxp   compact  form  based,  not  upon  the  mathematics  of  the  problem, 
but  upon  the  geometry  of  the  finite  element  model.   The  same  solution 
methods  may  be  used,  with  the  only  difference  being  that  the  inverting 
of   [A]   is  replaced  by  iteration.   This  iteration  can  only  be  as  accurate 
as  the  imposed  error  criterion.   For  large  systems,  however,  where  the 
[Aj    produces  a  full   NxN  matrix  on  the  RHS,  this  iteration  is 
especially  beneficial. 

The  last  contribution  was  the  development  of  a  new  equation-solver 
based  on  a  different  theory  than  the  DVOGER  routine.   The  new  code  deals 
with  the  differential  equations  directly  and  employs  the  Nxp   compacting. 


30 


VII.   TEST  PROBLEMS  AND  RESULTS 

A.  THE  PHYSICAL  MODEL 

A  cylindrical  nuclear  reactor  consisting  of  a  core  region  and 
blanket  having  overall  a  height  of  220  cm  and  a  diameter  of  180  cm  was 
modeled  for  the  test  problems.   As  shown  in  Figure  1,  a  radial  slice 
from  this  reactor  was  discretized  by  a  finite  element  of  54  elements  and 
38  nodal  points.   After  a  detailed  analysis  of  this  mesh  was  completed, 
preliminary  investigation  was  begun  using  a  model  of  220  elements 
(Figure  3). 

B.  COMPUTER  PROCESSING   CONSIDERATIONS 

All  programs  were  written  in  the  FORTRAN  IV  language  and  processed 
on  the  IBM  360/67  computer  using  the  FORTRAN  'H'  compiler.   The  results 
presented  in  terms  of  CPU  storage  requirements  and  processing  time  ought 
to  be  indicative,  in  a  relative  way,  of  those  obtained  from  other  machines. 
Single  precision  (seven  significant  digits)  was  used  throughout  the 
investigation,  with  the  results  later  verified  using  double  precision 
(fifteen  significant  digits).   The  observed  relative  discrepancy  was  less 
than  .01%.   It  is  recongized,  however,  that  for  larger  systems  a  more 
significant  difference  may  result. 

C.  PROBLEM  ANALYSIS 

The  three  techniques  discussed  have  been  incorporated  into  computer 
codes  and  have  been  tested  on  those  problems  described  in  Reference  3, 
namely,  the  dynamic  reactor  response  to   1)  a  disturbance  at  the  center, 

2)  a  uniform  disturbance  throughout  the  core,  and  3)  a  disturbance  occurring 

3 
at  a  skew  point  in  the  core   (R  =  40  cm,  Z  ■  0  cm).   The  direct,  or  N 
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2 
matrix  (an  Nx(N  )  array)  equation  was  solved  in  Reference  3,  and  those 

results  are  used  in  this  thesis  as  a  standard  for  comparison.   The 

computer  codes  developed  here  have  been  constructed  within  the  framework 

3 

and  parameters  of  the  original  N   program,  thereby  allowing  any 

discrepancy  to  be  attributed  to  the  technique  employed.   For  each  tech- 
nique and  disturbance,  the  flux  at  three  sample  points  was  recorded 
throughout  the  time  history  of  the  solutions  and  written  onto  a  computer 
storage  device.   A  separate  program  was  written  to  process  this  infor- 
mation and  produce  the  correlations  in  tabular  and  graphic  form.   The 
graphic  results  are  presented  in  Figures  4-21  for  the  center  point 
(R  =  0  cm,  Z  =  0  cm) ,  a  core  point  (R  =  40  cm,  Z  =  40  cm) ,  and  a  reflector 
point  (R  =  75  cm,  Z  =  80  cm).   Table  II  lists  the  specific  programs 
employed  arranged  by  disturbance  type   and  gives  parameters  such  as 
storage  requirements  and  computer  processing  time. 
1.    Solutions  with  the  CRANKO  Subroutine 

All  techniques  appear  to  give  acceptable  accuracy,  however 
there  is  a  startling  difference  in  processing  time.   The  CRANKO  subroutine 
was  able  to  accomplish  in  seconds  the  work  which  required  minutes  for  the 
DVOGER  subroutine  (See  Table  II).   The  steady  state  solution  given  by 
CRANKO  was  in  every  case  higher  than  that  obtained  by  DVOGER.   The 
implication  here  is  not  clear  since  there  is  no  analytical  solution 
available  to  represent  the  exact  answer.   In  general,  the  transient 
solutions  agreed  well,  except  for  the  oscillation  of  the  CRANKO  solution 
in  the  early  stages,  most  clearly  depicted  in  the  extreme  case  of  Figures 
8  and  9.   This  fluctuation  is  most  likely  the  result  of  inadequate  control 
of  the  CRANKO  time  step  selection  where  apparently  an  interval  has  been 
chosen  which  is  too  large  to  maintain  good  accuracy.   It  is  noteworthy, 
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however,  that  the  technique  is  so  stable  that  even  after  large  deviation 

it  is  able  to  recognize  the  error,  correct  itself,  and  return  to  the 

solution  curve  given  by  the  DVOGER  routine. 

2.    Solutions  with  the  DVOGER  Subroutine 

The  linearized  technique  using  premultiplication  by  the  previous 

time  step  flux  (\b      .  )   conformed  to  theoretical  expectations.   As  the 

t-At 

steady  state  solution  was  approached,  this  method  proceeded  with  in- 
creasing difficulty.   Hence  much  more  processing  time  to  advance  through 

the  steady  state  was  required. 

3 

The  direct   N   method,  on  the  other  hand,  has  no  difficulty 

in  recognizing  the  steady  state  and  marches  forward  with  large  time 

increments  until  problem  termination.   A  very  successful  effort  was  made 

to  use  the  exact  treatment  with  the  linearized  technique.   As  expected, 

more  processing  time  was  consumed  in  getting  to  the  steady  state  due  to 

the  requirement  to  recompute  the   [C  <p]   for  each  trial  ty,    but  once  the 

steady  state  was  approached,  the  procedure  progressed  rapidly  as  in  the 

direct  treatment.   Comparing  the  direct  with  the  linearized  exact  method, 

theory  developed  in  this  thesis  predicts  identical  results  with  twice 

the  computing  time  required  by  the  linearization.   As  shown  in  Table  II, 

this  is  indeed  the  case.   The  largest  relative  error  observed  between  the 

two  methods  was  .06%,  which  is  attributed  to  roundoff  steming  from  the 

doubling  of  the  computational  requirements.   When  the  greatest  accuracy 

is  required,  the  exact  treatment  with  linearization  may  therefore  prove 

the  most  useful,  especially  if  the  problem  rapidly  reaches  steady  state. 

This  discussion  of  the  linearized  technique  is  applicable  both  to  the  NxN 

method  or  to  the  Nxp   compact  method.   However  since  the  Nxp  method 

requires  iteration,  more  processing  time  and  less  accuracy  might  be 

expected. 
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D.    EVALUATION  OF  ERROR 

The  error  analysis  for  these  procedures  must  remain  subjective  as 
there  is  no  "correct"  result.   Several  sources  of  error  are  present  and 
deserve  comment.   Tests  were  made  with  both  EPS  =  .1  and  EPS  =  .01,  the 
error  criterion  for  the  DVOGER  routine,  and  little  difference  was  noted 

for  the  large  increase  in  processing  time.   [EPS  is  defined  by  DVOGER  as 

V   ERROR  (i)  ,    _.  r   .   .    ,   .   . 

I       — ~ * — *■  J.   The  convergence  criterion  [relative  deviation  of  any 

i=l     MAX 

\b.      between  successive  iterations]  for  the  iterative  solutions  in  CRANKO 
1 

was  varied  by  tenths  from  .01  to  .001  with  the  result  of  producing  higher 
steady  state  values  of  only  a  few  percent  as  the  criterion  was  made  more 
stringent.   Processing  time  for  this  variation  increased  only  a  few 
seconds.   Double  percision  trials  for  the  linearized  DVOGER  method  and 
the  compact  CRANKO  were  made  with  neglibible  difference. 

The  grid  of  only  54  elements  theoretically  introduces  the  largest 
source  of  error  from  the  "true"  solution.   With  the  techniques  developed 
in  this  thesis,  a  finer  grid  size  of  220  elements  has  been  implemented 
into  the  CRANKO  and  the  linearized  codes.   Because  the  linearized  code 
(DVOGER)  requires  about  566K  core  size,  extensive  testing  with  this 
technique  is  not  anticipated,  but  the  CRANKO  routine  which  requires  only 
188K  is  currently  being  investigated.   As  the  results  of  this  investiga- 
tion become  clear,  a  better  understanding  of  the  error  induced  by  each  of 
the  various  approaches  to  this  problem  would  be  obtained. 
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VIII.   CONCLUSION 

It  appears  that  nonlinear  differential  equations  resulting  from  a 
finite  element  formulation  (or  for  that  matter  finite  differencing)  may 
be  succesfully  reduced  using  the  linearizing  technique  developed  in 
this  thesis.   The  compacting  scheme  for  matrix  equations  effects  gross 
savings  for  storage  of  large  systems  (Table  I)  and  becomes  even  more 
attractive  when  the  inverse  of  a  matrix  is  either  not  desired  or  required. 
The  Crank-Nicolson  formulation  coupled  with  Gauss-Seidel  iteration  as 
expressed  in  the  computer  code  CRANKO  has  demonstrated  itself  to  be  a 
viable,  extremely  fast  technique  for  the  solution  of  differential 
equations.   Further  improvements  to  the  CRANKO  routine,  specifically  in 
regard  to  better  control  of  the  time  step  selection,  may  improve  upon  the 
attractiveness  of  this  approach. 

Development  is  warranted  in  several  areas.   The  computation  of  the 
Jacobian  using  iteration  must  be  formulated  in  order  to  use  the  DVOGER 
subroutine  with  the  compact  system.   DVOGER  itself  might  usefully  be 
reviewed  to  reduce  its  size  requirements  also.   Investigation  of  the 
problems  using  the  fine  mesh  size  of  220  elements  is  continuing.   This 
should  yield  more  exact  results  and  offer  a  better  standard  for  comparison. 
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APPENDIX  A 
TABLES 


TABLE   I 

IBM  360/67  CORE  REQUIREMENTS 
in  K  Bytes 


SYSTEM 

NO.  NODES 

OVERHEAD 

MATRIX  STORAGE 

TOTAL 

NxN2 

38 

100 

225 

325 

NxN 

38 

100 

35 

135 

Nxp 

38 

100 

12 

112 

NxN 

132 

150 

418 

568 

Nxp 

132 

150 

38 

188 

Nxp 

400 

290 

112 

402  t 

fcalculated   estimate 
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TABLE  II 

COMPARISON  OF  METHODS 
38  NODAL  PT.  MODEL 


METHOD 


PROCESSING  TIME  (MIN)      STEADY-STATE 
CORE  (k)    CENTER   UNIFORM   SKEW      ERROR  (%) 


DVOGER  N   .01 
DVOGER  N3  .1 
DVOGER  NxN  exact 
DVOGER  NxN  approx 
CRANKO 
CRANKO  DP 


330 
330 
130 
130 
110 
158 


5.32 

10.22 

10.73 

.62 

.77 


5.33 

6.27 

3.28 

4.57 

6.07 

8.10 

>15. 

>10 

.20 

.22 

.23 

.25 

<.01 
.06 
.16 

5. 

5. 


Notes 


1.  Processing  time  is  that  time  required  by  the  program  to  reach  one 
second  in  problem  time.   Steady  state  was  obtained  for  each  run. 

2.  .01  &  .1  pertain  to  DVOGER  error  criterion  EPS.   All  NxN  trials 
used  EPS  =  .1  CRANKO  used  iteration  convergence  tolerance  of  .0001. 

3.  Double  precision  (DP)  trials  showed  a  difference  in  the  fifth 
significant  digit  when  compared  to  the  SP  counterparts. 

3 

4.  NxN  approximate  technique  proceeds  as  rapidly  as  the  N  method 

until  steady  state  is  approached.   This  has  been  verified  by 
experiments  using  various  processing  time  cutoff  points. 
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APPENDIX  B 
FIGURES 
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REACTOR  MODEL  WITH  132  NODES 


t*Y\j 


/O  zo         30  40  SO         €0     65        -?&        85    &0 


FIGURE    3 
42 


GW 


y/o 


neutron 


'3     cm  -sec 


35" 


_.3o 


--  ZS 


ZO 


—  15 


-JO 


s 


-5 


loo    (sec) 


-3 


-2. 


-I 


TIME  DEPENDENT  NEUTRON  FLUX  AT  REACTOR  CENTER 

FOR  A  CENTRAL  DISTURBANCE 

PLOT  OF  CENTER  DISTURBANCE  -  CENTER  POINT 

FIGURE  4 


A3 
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DEVIATION  FROM  N  SOLUTION  AT  CORE  POINT  Z=40cm,  R=40cra 
COMPARISON  -  CENTER  DISTURBANCE  -  CORE  POINT 

FIGURE  7 
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DEVIATION  FROM  N   SOLUTION  AT  CORE  POINT  Z=40cm,  R=40cm 
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I.    DISCUSSION 

A.    The  FED- 2  Programs 

Both  of  the  FED-2  programs  consist  of  a  small  calling  program  which 
initiates  the  problem.   The  rest  of  the  program  is  stored  in  the  computer 
library.   The  advantages  of  a  pre-compiled  main  program  in  time  and 
convenience  is  clear,  but  the  use  of  dummy  dimension  statements  to 
pass  data  storage  throughout  the  calling  of  the  external  subroutine 
DVOGER  and  its  user-supplied  subroutine  appears  as  a  trick  not  normally 

considered  possible  with  the  FORTRAN  language. 

3 

When  the  original  N   FED-2  program  was  written,  every  available 

device  for  reducing  storage  was  attempted.   Thus,  where  properties  are 
normally  summed  over  areas  in  the  finite  element  method,  the  operations 
were  carried  out  on  the  nodal  points,  since  there  are  always  fewer  node 
points  than  elements.   The  new  FED-2  programs  continued  these  procedures 
not  only  to  provide  a  valid  comparison  but  also  because  this  economizing 
technique  is  felt  to  be  fundamentally  sound  even  if  unconventional.   Employ- 
ing this  device,  NFULEL  is  the  number  of  node  points  associated  with  the 
core  (fueled)  region,  whereas  it  would  conventionally  be  the  number  of 
core  (fueled)  elements. 

For  these  programs,  three  interest  points  IPTl,  IPT2,  and  IPT3,  may 
be  selected  for  detailed  investigation,  since  the  purview  of  information 
produced  by  the  printed  output  of  these  programs  is  quite  large,  making 
detailed  study  difficult.   The  time  history  of  the  neutron  flux  (PSI)  at 
each  of  these  interest  points  is  recorded  throughout  the  problem  in  an 
array  BATCH.   If  it  is  desired  to  use  this  collection,  it  can  be  written 
onto  a  computer  storage  by  the  inclusion  of  additional  JCL  in  the  GO  step. 
In  this  manner,  data  was  saved  for  analysis  by  the  data  processing  program 
used  in  this  thesis  to  prepare  the  graphs  of  Figures  4  -  21. 
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An  options  chart  and  a  schedule  for  FED-2  input  data  deck  listing 
are  included  in  this  section  in  addition  to  a  listing  of  the  programs. 
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OPTION  SELECTIONS  FOR  FED- 2  PROGRAMS 

NCOUNT 

=  0  PREMULTIPLICATION  BY  i> 

rt-At 

=  1  PREMULTIPLICATION  BY  $   TRIAL 

MTH 

=  1  DVOGER,  Jacobian  must  be  supplied 

=  2  DVOGER,  computes  its  own  Jacobian 

=  0  DVOGER,  doesn't  use  Jacobian 

=  5  CRANKO  routine  specified,  DVOGER  not  activated 


ERRVAL 


1-  used  internally  by  DVOGER,  card  input  value 
has  no  effect 

2-  used  by  CRANKO,  it  is  the  convergence  criterion 
for  satisfactory  ty   values 


EPSVAL 

1-  convergence  criterion  for  DVOGER 

2-  for  CRANKO,  sets  the  relative  problem  time 
difference  required  for  a  solution  to  be 
printed  out 
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c 
c 
c 
c 
c 

C         DATA  PROCESSOR  FOR  FED2  PROGRAMS. 
C 

c 
c 
c 

C    THIS  PROGRAM  TAKES  THE  STORED  DATA  SETS  AND  PRESENTS  A 

C    DIFFERENT  METHOD  AND/OR  DISTURBANCE.   THE  METHODS  MUST 

C    BE  IN  THE  ORDER  GIVEN.  ALL  METHODS  FOR  ONE  DISTURBANCE 

C    ARE  RUN  TOGETHER.   THE  CROER  FOR  DISTURBANCE  TYPES  IS 

C    1-CENTRAL,  2-UNFIROM,  3-SKEW.   THE  PROGRAM  IS  EASILY 

C    MODIFIED  TO  ACCEPT  CHANGES,  EG.  THE  EXCLUSION  OF  A 

C    DATA  SET  OR  METHOD, OR  THE  EXAMINATION  OF  A  DIFFERENT 

C    TIME  INTERVAL. 

C 

C  THE  METHODS  ARE 

C  EXACT      N**3 

C  LINEARIZED         N*N 

C  CRANK-NICQLSON 

C  GEAR         COMPACTED 

C  ADAMS      COMPACTED 

C 

C 

C         EXPLANATION  OF  ITEMS. 

C 

C    NIT  IS  THE  NUMBER  OF  TIME  POINTS  CCNTAINED  IN  A  GIVEN 

C    DSET.   EACH  TIME  A  DSET  IS  READ  IN,  ITS  CORRESPONDING 

C    NIT  MUST  BE  READ  IN  FROM  A  CARD. 

C 

C    JJ  IS  THE  NUMBER  OF  DIFFERENT  DISTURBANCES  EXAMINED. 

C 

C    JTYPE  IS  THE  NUMBER  OF  DIFFERENT  METHODS. 

C 

C  TCUT    IS    THE    LAST    NORMALIZED    TIME    POINT. 

C 

C    TSTRT  IS  THE  APPROXIMATE  FIRST  NORMALIZED  TIME  POINT. 

C 

C    NUM  IS  THE  NUMBER  OF  NORMALIZED  TIME  POINTS. 

C 

C 

C         HOUSKEEPING. 


C 


DIMENSION   FLUXC 200) ,FTIM(200) 

DIMENSION  STUF(50) ,T< 50 ), HOLD ( 400 ,4} , DSET (400,4) 
DIMENSION   EXACT (50,3) ,XLIN(50,3) ,ADAM(50,3), 
1    CRNK( 50,3) ,GEER(50,3) 
DIMENSION   DIFX( 50,3),DIFA(  50,3) ,DIFG( 50,3)  ,DIFC(50,3) 

DIMENSION  P(50) 

DIMENSION  H( 1800) 

OIMENSION  STOF.(50)  ,ITB(12),RTB(28) 

DATA  ITB,RT3/12*0,28*.0/ 

DATA  CTRrUNF ,SKW/«CENR« , 'UNIF' , 'SKEW1/ 

DATA  CTRL, CORE, REEL/' CTRL' , « CORE • , 'REEL • / 

OATA  CR/' CRNK' / 

DATA  GE/'GEAR'/ 

DATA  AD/'ADAMV 

DATA  EX/'EXCT'/ 

DATA  XL/'XLIN'/ 

DATA  PT/«  ....'/ 

110  FORMAT  (5X,I5) 

205  FCRMAT  (•  1  •,  10X  ,'  S  INGLE  CENTRAL  DISTURBANCE'/) 

207  FORMAT  ( • 1 ' , 10X , ■ UNI F ORM  DISTURBANCE'//) 

209  FORMAT  ( •  1 '  , 10X , ' SKEW  DISTURBANCE1//) 

211  FORMAT  (10X, ' EXACT'//  ) 

212  FORMAT  (10X, 'XLIN' //) 

213  FORMAT  (10X, 
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c 
c 

c 
c 
c 


c 
c 
c 

c 


c 
c 
c 
c 
c 
c 


c 
c 
c 


214 
215 
222 


224 
232 

234 
236 
238 
240 
242 

i 
244 
246 

1 
5CC 


04 
03 


21 


FCRMAT 
FORMAT 
FORMAT ( 

•PSI 

•TIM 
FORMAT 
FORMAT 
FCRMAT 
FORMAT 
FCRMAT 
FCRMAT 
FORMAT 

•GEE 
FORMAT 
FORMAT 

8X,F 
FORMAT 


ESTABLISH    THE    CONTROLING    PARAMETERS. 

SCALE    =1.E    14 

TSTRT=1 .E-5 

NUM=50 

VAL=.05 

JTYPE=3 

JJ  =  3 


ZERO    THE    STORAGE    AREAS. 


DO    03 

DO    04 

EXACT 

XLIN 

A  CAM 

GEER 

CRNK 

DIFX 

DIFA 

DIFG 

DIFC 

CONTINUE 

CONTINUE 


1=1 
J  =  l 
(  It 
(  I, 
(  I, 
(  It 
(  It 
(  It 
(  I, 
1  It 
(  It 


NUM 

5 

)  =  .0 

)  =  .0 
)  =  .0 
)=.0 
)  =  .0 
)=.0 
)  =  .0 
)  =  .0 
)  =  .0 


DC    88    LL=1,JJ 

LL    DENOTES    WHICH 
CONSIDERED. 


TYPE    OF    DISTURBANCE     IS    BEING 


CONSTRUCT    THE    NORMALIZED    TIME     INCREMENTS 

XNUM=NUM 

TCUT=.l 

IF    (LL     .EG.    2)       TCUT=.0329 

X=(TCUT/TSTRT)**( 1 ./XNUM) 

DO    21    1=1, NUM 

T(I )=TSTRT*(X**I  ) 

P(I)=    ALCG101TC I  )  ) 

CONTINUE 

DO    86    LM=1 ,JTYPE 

REAC    AND    WRITE    ONE    DATA    SET     IN    THE    RAW. 


READ    (5, 110)    NI 
READ    (4)    DSET 
IF    (LL     .EC.     1) 
IF    (LL     .EQ.    2) 
IF    (LL     .EQ.    3) 
IF    (LM  .EQ. 


WRITE(6,205) 
wRITE(6t207) 
WRITE(6,209) 
1)     WRITE(6,2il) 
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c 
c 
c 
c 
c 


c 

c 

c 
c 


c 

c 


c 
c 
c 
c 
c 


c 
c 
c 
c 
c 


c 
c 
c 


06 


05 
07 


23 


25 


26 


13 
15 

19 


30 


IF 
IF 
IF 
IF 
IF 


(LM 
(LM 
(LM 
(LM 
(LA 


WRITE  (6 


2) 
3) 
3) 
4) 
5) 


WRITE(6,212) 
WRITE(6,213) 
WRITE(6,215) 
WRITE(6t2l4) 
WRITE(6t 215) 


.EQ. 
.EQ. 
.EQ. 
.EQ. 
.EQ. 
224)     ( (DSET( I ,K),K=1,4) ,1=1, NIT) 


REBUILD    THE    DATA 
ARE     IDENTICAL. 


SET    EXCLUDING    TIME    VALUES    WHICH 


DC    06    1=1,4 

HOLD(  1,  I)=    DSETd, 

CONTINUE 


I) 


JK=1 

PC    07    1=2, NIT 

J  =  JK 

IF    ((DSET( I,1)-H0LD( J,l) )/HCLD( J,l)     .LT.     VAL)     GO    TO    07 

JK=JK*1 

DO    05    JIM=1,4 

HOLD( JK, JIM)=DSET( I , JIM) 

CCNTINUE 

CONTINUE 


DO    17    M=2,4 


EXCLUDE 
rJK 


*    1 


SCREEN    TO 
DO    23    JEND=1 
JDIF=JK-JEi\D 
JP  =  JOIF-H 
IF    (HOLD(JDIF, 
CONTINUE 
JP=JDIF 
GO    TO    26 
CONTINUE 
DO    26     1=1, JP 

SMOOTH    OUT 
HOLD( I ,M)=HQLD( JP, 
CCNTINUE 


DATA    POINTS    GOING    NEGATIVE 


^)     .LT.     .0)     GO    TO    25 


ANY 


DISCONTINUITY. 

M) 


DO  27  1=1, NUM 
STUF(  I)=P(  I) 
27  CCNTINUE 


INSURE  THE  LAST  TIME  VALUE  FOR 
LAST  NORMALISED  TIME  POINT. 


DSET  IS  BEYOND  THE 


JB  =  0 

DO  15   JF=1,JK 

JA=JK-JB 

JB=JB+1 

IF  <H0LD(JA,1)  .GE.  T(NUM))   GO  TO  19 

DO  13  1=2,4 

HOLDUA,  I)=SCALE 

CCNTINUE 

CCNTINUE 

CCNTINUE 


SCALE  THE  INPUT  VALUES  FOR  THE  INTERPOLATION  ROUTIN 

DO  30  1=  1,JK 

FT IM< I ) =AL0G10(H0LD( 1,1)) 

FLUX(I)=       (HOLD( I,M) )/SCALE 

CONTINUE 

WRITE  (6,224)  ( ( HOLD( I ,K) , K=l , 4 ) , 1=1 , JK ) 


109 


c 
c 


c 
c 
c 
c 
c 

c 

c 


NORMALIZE  THE  DATA  SET 


C 
C 
C 
C 
C 


41 


42 
43 


44 
45 


46 
47 


48 
49 


50 
51 
17 
86 


61 
63 


JER=0 
JT  =  JK 
CALL    ICS1VU 


(FLUX,FTIM,JT,NUM,STUF,H, JER) 


TRANSFER  NORMALIZED  DSET  TO    APPROPRIATE  CLASS  AND 
COMPUTE  THE  DIFFERENCE. 

L=M-1 

L  ESTABLISHES  WHICH  POINT  IN  THE  REACTOR  IS  BEING 

EXAMINED. 
GO  TO  (41,  43,  45,  47,  49), LM 
CONTINUE 
DO  42  1=1, NUM 
EXACT(I ,L)=STUF( I) 
CONTINUE 
GO  TO  51 
CONTINUE 
DG  44  1=1 
DIFX(I,L) 
XLINd  ,L)=STUF(  I  J 
CONTINUE 
GO  TO  5  1 
CONTINUE 
DO  46  I=1,NUM 
CPNK( I, L)=STUF( I) 
DIFC( I,L)=((STUF( I ) 
CONTINUE 
GO  TO  51 
CONTINUE 
DO  48  1=1, NUM 
DIFG( I, L)=((STUF( I ) 
GEER( I,L)=STUF( I ) 
CCNTINUE 
GO  TO  51 
CCNTINUE 
DC  50  1=1, NUM 
DIFA(  I,L)  =  ((  STUF( I ) 
ACAM( I,L)=STUF( I) 
CCNTINUE 
CCNTINUE 
CCNTINUE 
CCNTINUE 

DO  63  L=l,3 

DO  61  1=1, NUM 

EXACT ( I ,L)=EXACT( I ,L)*SCALE 

CCNTINUE 

CONTINUE 


,  NUM 
(( STUF(  I  )  -EXACT ( I,L) J/EXACTC I ,LJ  )*100. 

-SCALE 


*SCALE 
-EX ACT (  I,L) ) /EXACT ( I ,L) )*100. 


-EX ACT (  I,L)  ) /EXACT (I  ,L)  )*100. 
*SCALE 


-EXACT ( I,L) ) /EXACT (I ,L)  )*100. 
*SCALE 


DO 
DO 
IF 
IF 
IF 


WRITE    OUT    RESULTS    AND    DRAW    THE    GRAPHS 
57    LCNT=1,4 


58  L    =1,3 

(LL  .EQ.    1) 

(LL  .EC.    2) 

(LL  .EQ.    3) 


WRITE(6,205) 
WRITE(6,207) 
WRITE( 6,209) 


WRITE  (6,240) 
IF(L  .EQ.  1) 
IF(L  .EQ.  2) 
IF(L  .EQ.  3) 
IF    (LCNT     .EQ. 


WRITE(6,232) 
WRITE(6,234) 
WKITE(6,236) 
2)    WRITE16,238) 


WPITE(6,242) 

IF    <LCNT     .EQ.l)     WRITE(6,244) 

1  (T(I ) , EXACT ( I,L),XLIN(I,L),GEER(I,L),ADAM(I,L) ,CRNK 

2  1=1, NUM) 
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c 
c 
c 


IFUCNT     .E 

1  (T(I )  ,E 

2  1=1, NUM 
LPLT=LPLT+ 
IF  (LCNT  . 


Q.  2)  WRITE(6,246) 

XACTU  ,L),DIFX(I,L),DIFG(I,L),DIFA(I,L),DIFC 

1 

LE.  2)  GO  TO  54 


INSERT  CALLS  TO  PLOTTING  ROUTINES. 


C 

c 


71 


72 
73 


74 
75 


76 
77 

78 
79 


80 
95 

97 

54 
5c 
57 


RT8(6)=PT 
RTB(2)=.0 
IF  (LCNT 
(LL  . 
(  LL  . 
(LL  . 
(L 
(L 
(L 


IF 
IF 
IF 
IF 
IF 
IF 


EQ. 
EQ. 
EQ. 


DC  97 

GO  TO 

CCNTI 

RTB(3 

ITB(1 

DO  72 

STUF( 

IF  (L 

CONTI 

GO  TO 

CONTI 

RTB(3 

ITB(  1 

DC  74 

STUF( 

IF  (L 

CCNTI 

GO  TO 

CCNTI 

RTB(3 

ITBU 

DC  76 

STUFi 

IF  (L 

CONTI 

GC  TO 

CONTI 

RTB(3 

DO  78 

STUF( 

IF  (L 

CCNTI 

GC  TO 

CCNTI 

RTS(3 

DC  80 

STOF< 

IF  (L 

CCNTI 

CCNTI 

CALL 

CGNTI 


EQ.  4)  RTB(2)=5. 
Q.  1)   RTB(5)   =  CTR 
C.  2)   RTB(5)   =UNF 
Q.  3)   RTB(5)  =  SKW 

1)  RTB(7)    =CTRL 

2)  RTB(7)    =CORE 
3J  RTB(7)    =REFL 

T=l, JTYPE 

73, 75,77,79), IPLOT 


NUM 

EXACT( I ,L) 

EQ.  4)  STUF( I)=0. 


NUM 

I N ( I , L) 

EQ.  4)  STUF(  I)  =  DIFX(  I,LJ 


=  CR 


I  PLO 

(71, 
NUE 
)   =E 
)  =1 

I  =1  • 
I)  = 
CNT  . 
NUE 

95 
NUE 
)   =X 
)  =  2 

1=1, 
li  =XL 
CNT  . 
NUE 

95 
NUE 
) 
)  =  3 

1  =  1, 

n  =cr 

CNT  . 
NUE 

95 
NUE 
I   = 

1  =  1 
I)=GE 
CNT  . 
NUE 

95 
NUE 
)   = 

I=l» 
I)  = 
CNT  . 
NUE 
NUE 

DRAWP(NUM,P,STUF,ITB,RTB) 
NUE 


NUM 

NK(I,L) 

EQ.  4)  STUF( I)=DIFC( I,L) 


GE 

,  NUM 

ER(I ,L) 

EQ.  4)  STUFi I)=DIFG( I, L) 


AD 

NUM 

ADAM( I ,L) 

EQ.  4)  STUFt I)=DIFA( I,L) 


CCNTINUE 
CCNTINUE 
CCNTINUE 


88  CGNTINUE 

WRITE  (6,5 

STOP 

END 


00) 
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c 
c 

C  A  PROGRAM  TG  TEST  THE  COMPACT  METHOD 

C  OF  FINITE  ELEMENT  MATRIX  WORK. 

C 

C 

C 

C 

C    INPUT  INFORMATION 

C 

C    DATA  CARDS  ARE  REQUIRED  IN  THE  FOLLOWING  ORDER, 

C 

C    NODY,  NUMSNP   ( NR  NODES  IN  Y  DIRECTION,  NR  SYSTEM 

C   NCDAL  POINTS) 

C 

C    PSI   THE  SET  OF    VALUES  FOR  TEST  MULTIPLICATION. 

C 

C    NUMEL    NR  ELEMENTS  IN  THE  MESH 

C 

C    ELNOD    THE  ELEMENT  CONNECTIVITY  MATRIX. 

C 

C 

c 

INTEGER*^  ELNOD 

DIMENSION  VALUE ( 50) ,MN0D(50,7) , VMATX(50, 7) , 
i    BIGA( 50,50) ,ELNGD( 75,3) 
DIMENSION   PSK50),  COMPCT( 50 ) , USUAL( 50 ) 
DIMENSION  AA(3,3) 
110  FORMAT  (I15,5X,I5) 

119  FORMAT  (« l«/////////«      •) 

120  FORMAT  ( F9 . 2 , 6F1 0 . 2 ) 
150  FORMAT  (5X,3U0) 

C 

210  FORMAT*  ■  J.1  ,10X,'N0DE  NEIGHBOR  CONNECTIVITY'//) 

212  FORMAT  ( 10X, • NODE • , 5X , • K2* , 5X, • K3 ■ , 5X , « K4« , 5X , • K5 ' , 

1    5X,« K6» ,5X, 'K7' //) 
215  FORMAT  ( 1  OX, 7(  I  5 , 2X) / ) 

220    FCRMATt • 1 'Z//10X, • THIS     IS    THE    NUMSNP    X    7    MATRIX'//) 
225    FORMAT     ( 1  OX, 7( F6 . 0 , IX ) / ) 

230    FORMAT     ( •  1  •  ///10X , *  THE    ELEMENT    CONNECTIVITY    MATRIX1/) 
235    FORMAT     ( 1  OX , I  4, 3 ( 3X , I  4 ) ) 
250    FORMAT* ' 1 »///10X , ' THIS    IS    THE    USUAL    NUMSNP    X    NUMSNP    • 

1  ■BIGA    MATRIX' // ) 

2  55    FORMAT     (2X  ,  20 ( F5 . 0  ,  IX ) /) 

270    FORMAT     (• 1 «//// 1  OX , c COMPACT    FORM • ,5X, ' USUAL    FORM'//) 
272    FORMAT     ( 5X , I  2 , 7X , F8. 0, 07X , F8 .0 ) 
C 
C 

READ(5,110)     NODY, NUMSNP 
READ    (5,120)     (PSI        ( I ) , 1=1 , NUMSN P) 
WRITE    ( 6,119) 

WRITE    (6,120)     (PSI       (I),     1=1, NUMSNP) 
C 

C         CONSTRUCT    THE    NEIGHBOR    NODES    OF    THE    HEXAGON. 
WRITE    (6,210) 
WRITE     (6,212) 
C 

DO    35    Kl  =  l, NUMSNP 
N  =  l 
5    CONTINUE 
NN=N*NODY 
N  =  N*l 

IF    (( Kl-NN)     .GT.     0)       GO    TO    5 
K3=K1-1 
K6  =  KH-1 
K2=-Ki-NCDY 
K7=K2+1 
K5=K1+N0DY 
K4=K5-1 

IF(K1     .GT.    NODY)       GO    TO   6 
K2=0 
K7  =  0 
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6  IF(K1     .NE.    NN)       GO    TO    7 
K6=0 

K7  =  0 

7  IF(K1     .LE.     (NUMSNP-NODY))       GO    TO    8 
K5=0 

K4  =  0 

8  IF    (Kl     .NE.     (NN-NODY+1))       GO    TO    9 
K3=0 

K4  =  0 

9  CONTINUE 
C 

MNODCK1 ,  1)=K1 

MNOD(Kl ,2)=K2 

MNOD(Kl ,3)=K3 

MNOD(Kl ,4)=K4 

MNODtKl ,5)=K5 

MNGDiKl ,6)=K6 

MNOD(Kl ,7>=K7 

WRITE    (6,215)       (MNOD(Kl,I) ,1=1,7) 
35    CONTINUE 
C 
C 
C 
C  READ    THE    CONNECTIVITY    MATRIX 

READ    (5,110)    NUMEL 

WRITE    ( 6,230) 

DO    60       1=1 ,NUMEL 

READ    (5,150)     (ELNODd  ,  J)  ,  J=l,3) 

WRITE     (6,235)        I , ( ELNOD( I , J ) , J=l ,3) 

60  CONTINUE 
C 

C  ZERO    THE    3IGA,     COMPCT,     AND    VALUE    MATRICIES. 

DO   62    1=1 , NUMSNP 
VALUE( I )=0. 
DC    61    J=l ,NUMSNP 
VMATXU ,  J)=0. 
BIGA( I , J)=0. 

61  CONTINUE 
6  2    CONTINUE 

C 

c 
c 
c 

C      CCNSTRUCT  THE  N  X  7  MATRIX 
C 

DO  16  L=1,NUMEL 

Y=FLQAT(L) 

Ri=2.*Y 

R2=3.*Y 

R3=4.*Y 

AA(1, 1) =2.*Ri+R2+2.*R3 

AA( 1,2)=RH-R2+2.*R3 

AAtlt3)=2.*Rl+2.*R2+R3 

AA(2, i)=2.*RH-2.*R2+2.*R3 

AA(2,2)=R1+R2+R3 

AA(2,3)  =  RH-2.*R2  +  2.*R3 

AA(3, 1)=AA( 1,1) 

AA(3,3)=AA(2,3) 

AA(3,2)=AA(2,1) 

DO  15  K=l,3 

KK=ELNOD(L,K) 

DO  14  1=1 ,1 

II=ELNOD(L,I ) 
C 

C      SEARCH  KK  ROW  FOR  NODE  II.   IF  FOUND,  POST  THE  VALUE 
C      TO  VMATX(KK,M).   OTHERWISE,  CONTINUE  STEPPING 
C      THROUGH  THE  ELEMENT. 

DO  12  M=l,7 

NOW  =  M 

MM    =MNOD(KK,M) 

IF  (MM     .EO.  II )  GO  TO  13 
12  CONTINUE 
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13  CGNTINUE 
VMATX(KK,NGW)=VMATX(KK,NOW)  +  AA(K,I) 

14  CGNTINUE 

15  CGNTINUE 

16  CCNTINUE 
C 

C      WRITE  OUT  THE  N  X  7 

WRITE  (6,220) 

WRITE  (6,212) 

DO  37  K1=1,NUMSNP 

WRITE  (6,225)  ( VMATX ( Kl , J ) , J  =  l , 7 ) 
37  CCNTINUE 
C 
C 
C 

C      GENERATE  THE  ELEMENT  VALUES  AND  MOVE  THEM  INTO  THE  NXN 
C 

DO  68  L  =  l ,NUMEL 

Y=FLOAT(L) 

R1=2.*Y 

R2=3.*Y 

R3=4.*Y   • 

AA(1,1)=2.*RH-R2  +  2.*R3 

AA( 1,2) =R1+R2+2.*R3 

AA(  1,3)=2.*RH-2.*R2  +  R3 

AA(2,1 ) =2.*Rl+2.*R2+2.*R3 

AA(2,2) =R1+R2*R3 

AA(2,3) =R1+2.*R2+2.*R3 

AA(3,1) =AA( 1,1) 

AA(3,2)=AA(2,1 ) 

AA(3,3) =AA(2,3) 

DO  67  K-1,3 

KK=  ELNOD(L,K) 

DC  66  1=1,3 

I I=ELNOD(L,I  ) 

BIGA(KK,II)=BI3A(KK,II)«-AA(K,n 

66  CONTINUE 

67  CCNTINUE 

68  CCNTINUE 
C 

C      WRITE  OUT  THE  N  X  N. 

WRITE  (6,250) 

DO  75  I=1,NUMSNP 

WRITE  (6,255)  (BIGA(I,J),  J=1,NUMSNP) 
75  CONTINUE 
C 
C 

c 

C      PERFORM  A  SIMPLE  MATRIX  X  COLM  MULTIPLICATION  ON  BOTH 

C      AND  COMPARE  THE  RESULTS. 

C 

DO  72  1=1 ,NUMSNP 

CCMPCT( I)=0. 

DO  71  J =1,7 

NAME=MNOD( I, J) 
COMPCT(I)=PSI(NAME)*VMATX(I,J) *  COMPCTt I) 

71  CCNTINUE 

72  CCNTINUE 
C 

DO    42    I=1,NUMSNP 

USUAL( I )=0. 

DO    41    J=l ,NUMSNP 

USUAL ( I  )  =  8IGA(  I, J)*PSI(J)+USUAL(  I  ) 

41  CONTINUE 

42  CONTINUE 


C 

c 


WRITE(6,270) 

DC    80    I=1,NUMSNP 

WRITE(6,272)     I,     COMPCT ( I ) , USU AL ( I ) 

80    CONTINUE 
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STOP 

END 

//GO, 

.SYSIN  DO 
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21 
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D. 
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c 

CMPACT  FORM 

USUAL  FORM 

1 

780. 

780. 

2 

8130. 

8130. 

3 

21120. 

21120. 

4 

25810. 

25810. 

5 

15799. 

15799. 

6 

54965. 

54965. 

7 

98365. 

98365. 

8 

72195. 

72195. 

9 

39740. 

39740. 

10 

126275. 

126275. 

11 

195386. 

195386. 

12 

113725. 

113725. 

13 

92605. 

92605. 

14 

24155o. 

241556. 

15 

288460. 

288460. 

16 

146967. 

146967. 

17 

74300. 

74300. 

18 

144144. 

144144. 

19 

144310. 

144310. 

20 

39144. 

39144. 
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THIS  IS  THE  NUMSNP  X  7  MATRIX 


NODE 


K2 


KB 


K4 


K5 


K6 


K7 


15. 

0. 

0. 

0. 

13. 

14. 

0. 

96. 

0. 

15. 

48. 

81. 

54. 

0. 

171. 

0. 

48. 

114. 

147. 

90. 

0. 

141. 

0. 

80. 

180. 

108. 

0. 

0. 

151. 

18. 

0. 

0. 

105. 

152. 

44. 

458. 

74. 

148. 

246. 

279. 

214. 

95. 

618. 

134. 

216. 

312. 

345. 

276. 

149. 

345. 

96. 

284. 

378. 

216. 

0. 

0. 

433. 

98. 

0. 

0. 

195. 

338. 

203. 

93  8. 

254. 

352. 

44<t. 

477. 

400. 

257. 

1098. 

314. 

420. 

510. 

543. 

462. 

311. 

549. 

192. 

488. 

576. 

324. 

0. 

0. 

709. 

182. 

0. 

0. 

285. 

524. 

36  5. 

1418. 

434. 

556. 

64  2. 

675. 

586. 

419. 
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